NDVI-HGM-WQ Model

Data Filtering and Merging

Take original data from previous and munge and merge with updated Dave Weixelman meadows that were selected.

wx.mdw<-read_rds("~/Documents/github/meadow_hgm_waterchem/data_output/wx_mdw_dat.rds") # weixelman data for 114 meadows

# rename and sort columns a bit
# Drop HGM source_type 8 (Dry), and combine classes 6 and 7, add elev col

wx.mdw <- wx.mdw %>% 
  select(PLOT, UCDavisObject_ID:PLOTNAME2, AREA_ACRE:Shape_Area) %>% 
  rename(Method=METHOD.y, UCDID=UCDavisObject_ID) %>% 
  mutate(hgm_class_comb = ifelse(source_type==8, NA, source_type), 
         hgm_class_comb = ifelse(source_type==7, 6, hgm_class_comb),
         elev_av_2100 = ifelse(ELEV_MEAN<2100, "<2100",">=2100"))

## make a quick meadows type table
mdw_types<-c("lacustrine fringe", "depressional",
             "discharge-slope-hillslope",
             "riparian-discharge-slope", "riparian",
             "subsurface-discharge-slope",
             "subsurface", "dry")
mdw_type_code<-c("lf", "dep", "ds", "rip-ds",
                 "rip", "sub-ds", "sub", "dry")
mdw_type_class<-seq(1,8,1)

# bind df
mdw_hgms<-data.frame("hgm_type"=mdw_types,"hgm_code"=mdw_type_code, "hgm_class"=mdw_type_class )

# now join this with original dataset
wx.mdw2<-inner_join(wx.mdw, mdw_hgms, by=c("source_type"="hgm_class"))
head(as.data.frame(wx.mdw2[,c(2,14,16,44:47)]))

# save for easy reload next time:
write_rds(x = wx.mdw2, path ="~/Documents/github/meadow_hgm_waterchem/data_output/wx_mdw_dat2.rds")

First filter data to only selected data (OK and not DRY or NOT OK) based on Dave Weixelman’s excel file from June 05, 2016. Several meadows were not accessible and were dropped, final list only includes data where there were UCDSNM_ID matched. Second, recombine and drop HGM classes (drop dry category 8 and merge 6 & 7). Finally write back out as .rds file.

wx.mdw<-read_rds("~/Documents/github/meadow_hgm_waterchem/data_output/wx_mdw_dat2.rds")

# easy ggplot
# wx.mdw %>% ggplot + aes(x=hgm_class_comb, fill=elev_av_2100) + geom_bar()

# this is the updated data version: 
wx.mdw.sub<-read_csv("~/Documents/github/meadow_hgm_waterchem/data/water_source_plot_selection_2016_06_05.csv")
head(wx.mdw.sub)

# filter out the "NOT OK" and select cols of interest:
wx.mdw.sub <- wx.mdw.sub %>% 
  mutate(selected=ifelse(Notes=="NOT OK" | Notes=="DRY", 0, 1)) %>% 
  filter(selected==1) 
wx.mdw.sub %>% tally(selected) # should be 26 here!

# keep only selected and see how many matches occur:
dim(inner_join(wx.mdw, wx.mdw.sub, by="PLOT"))  # only 19 match up!
dim(inner_join(wx.mdw, wx.mdw.sub, by=c("UCDID"="UCDavisObject_ID")))  # 24 match up!

# this is the refiltered data...based on Dave's comments. Re-merge with dataset
wx.mdw.sub.UCID<-inner_join(wx.mdw, wx.mdw.sub[,c(3,5,11,13)],
                            by=c("UCDID"="UCDavisObject_ID"))
wx.mdw.sub.UCID[2,48]<-6 # just make this a 6

# view the categories, note there are a few mismatches...ask Dave?
wx.mdw.sub.UCID %>% select(source_type.x,source_type.y, hgm_type, hgm_classes) %>% 
  as.data.frame()

# now re-format the hgm_class_comb column
wx.mdw.sub.UCID <- wx.mdw.sub.UCID %>% 
  mutate(hgm_class_comb = ifelse(source_type.y==8, NA, source_type.y), 
         hgm_class_comb = ifelse(source_type.y==7, 6, hgm_class_comb))

# ok write this out to merge with GEE data later
write_rds(x = wx.mdw.sub.UCID, path ="~/Documents/github/meadow_hgm_waterchem/data_output/wx_mdw_sub_UCDID.rds")

Merge with Google Earth Engine Data

Need to merge with original GEE data after removing duplicated columns, then filter, rename a col, and save out.

# HGM data
wx.mdw.sub<-read_rds(path = "~/Documents/github/meadow_hgm_waterchem/data_output/wx_mdw_sub_UCDID.rds")

# GEE aata
load(file = "data/mdw_geedat2.rda")

# Make a list of duplicated column names to filter from merge:
wx<-names(wx.mdw.sub)
gee<-names(mdws)
allcols<-c(wx, gee)
coldups <- allcols[duplicated(allcols)]

# now select columns that aren't duplicated before merging
library(magrittr) # to allow for special pipes %<>%
wx.mdw.sub %<>% select(-one_of(coldups), -source_type.x) # writes back to same object

# merge with GEE data
mdw.mod.dat<-inner_join(wx.mdw.sub, mdws, by=c("UCDID"="ID"))
mdw.mod.dat %<>% rename(source_type = source_type.y)
names(mdw.mod.dat)

write_rds(mdw.mod.dat, "~/Documents/github/meadow_hgm_waterchem/data_output/mdw_mod_GEE.rds")

Final Dataset

The Weixelman dataset originally selected 24 meadows, of which there were only 17 unique UCD Meadow IDs, as several of the meadows (PLOTNAME2) actually fall within the same UCD ID. In addition, the original Google Earth Engine (GEE) dataset that we used did not contain a number of the UCDID’s, probably because that data wasn’t able to be pulled via the GEE. Need to follow up with Andy on that.

# HGM data
wx.mdw.sub<-read_rds(path = "~/Documents/github/meadow_hgm_waterchem/data_output/wx_mdw_sub_UCDID.rds")

unique(wx.mdw.sub$UCDID)

##  [1] "UCDSNM012349" "UCDSNM012863" "UCDSNM005339" "UCDSNM000426"
##  [5] "UCDSNM001115" "UCDSNM000195" "UCDSNM016565" "UCDSNM016589"
##  [9] "UCDSNM015782" "UCDSNM012985" "UCDSNM013064" "UCDSNM016769"
## [13] "UCDSNM011600" "UCDSNM009725" "UCDSNM011455" "UCDSNM011447"
## [17] "UCDSNM014209"

unique(wx.mdw.sub$PLOTNAME2)

##  [1] "ELD0101 -ONION VALLEY"         "ELD0106 -GOVERNMENT MEADOW"   
##  [3] "INY0105 -UPPER BAKER"          "INY0113 -Bakeoven Meadow"     
##  [5] "INY0114 -Upper Monache Meadow" "INY0115 -Lower Monache Meadow"
##  [7] "INY0130 -POISON"               "INY0131 -ROUND"               
##  [9] "INY0132 -ADMIN SITE"           "INY0137 -MULKEY SOUTH"        
## [11] "LAS0103 -Upper Long Valley"    "LAS0104 -Lower Long Valley"   
## [13] "LAS0113 -Pasture 172"          "LTB0102 -UPPER ANTONE"        
## [15] "LTB0103 -ANTONE"               "MOD0103 -Cottonwood Spring"   
## [17] "SEQ0114 -Beck Meadow"          "SEQ0115 -Snake Creek Meadow"  
## [19] "STA0102 -KENNEDY LAKE"         "STA0103 -NIGHTCAP"            
## [21] "STA0104 -Indian Springs"       "STA0108 -Upper Relief Valley" 
## [23] "STA0109 -Hay Meadow"           "TAH0123 -BARKER MEADOW WEST"

See below for stats. The final GEE dataset merged with Weixelman selected meadows gives us 9 meadows with unique UCDIDs. Fine for first cut of model. In comparing NDVI and NDWI, it appears only 5 HGM groups were retained, the subsurface HGM type didn’t have any meadows selected and thus doesn’t exist in this datase. See table for description of types from original through final GEE dataset.

library(DT)

mdw_types<-c("lacustrine fringe", "depressional",
             "discharge-slope-hillslope",
             "riparian-discharge-slope", "riparian",
             "subsurface-discharge-slope",
             "subsurface", "dry")
mdw_type_code<-c("lf", "dep", "ds", "rip-ds",
                 "rip", "sub-ds", "sub", "dry")
mdw_type_class<-seq(1,8,1)

# bind df
mdw_hgms<-data.frame("hgm_type"=mdw_types,"hgm_code"=mdw_type_code, "hgm_class"=mdw_type_class )

DT::datatable(mdw_hgms, rownames = FALSE, autoHideNavigation = T,
              class="compact row-border stripe",
              caption=htmltools::tags$caption(
                style = 'caption-side: bottom; text-align: center;',
                htmltools::em('Table 1. '), 
                htmltools::em('All HGM Types used from Weixelman et al. 2011')
              ), colnames = c("HGM Type"=1, "HGM Code"=2, "HGM Class"=3)
) %>% DT::formatStyle('HGM Type',  color = 'black', 
                  backgroundColor = 'lightblue', fontWeight = 'bold')

mdw_hgm_rev<-mdw_hgms[c(1:5,7),]
mdw_hgm_rev[6,3]<-6
mdw_hgm_rev[6,3]<-6

DT::datatable(mdw_hgm_rev,rownames = FALSE,
              autoHideNavigation = T,
              class="compact row-border stripe",
              caption=htmltools::tags$caption(
                style = 'caption-side: bottom; text-align: center;',
                htmltools::em('Table 2. '), 
                htmltools::em('Revised HGM Types for selection of meadows')
              ), colnames = c("HGM Type"=1, "HGM Code"=2, "HGM Class"=3)
) %>% DT::formatStyle(columns = 'HGM Type',  color = 'black',fontWeight='bold') %>% 
  DT::formatStyle(columns = 'HGM Class',color = 'white', fontWeight = 'bold',
                  backgroundColor = styleInterval(cuts=c(1,2,3,4,5),
                                                  values = c("#440154FF", "#414487FF",
                                                             "#2A788EFF", "#22A884FF",
                                                             "#7AD151FF", "#FDE725FF")))

mdw_hgm_rev2<-mdw_hgm_rev[c(1:5),]
DT::datatable(mdw_hgm_rev2,rownames = FALSE,
              autoHideNavigation = T,
              class="compact row-border stripe",
              caption=htmltools::tags$caption(
                style = 'caption-side: bottom; text-align: center;',
                htmltools::em('Table 3. '), 
                htmltools::em('Revised HGM Types in final meadow dataset')
              ), colnames = c("HGM Type"=1, "HGM Code"=2, "HGM Class"=3)
) %>% DT::formatStyle(columns = 'HGM Type',  color = 'black',fontWeight='bold') %>% 
  DT::formatStyle(columns = 'HGM Class',color = 'white', fontWeight = 'bold',
                  backgroundColor = styleInterval(cuts=c(1,2,3,4),
                                                  values = c("#440154FF", "#3B528BFF",
                                                             "#21908CFF", "#5DC863FF",
                                                             "#FDE725FF")))

mdw.dat<-read_rds("~/Documents/github/meadow_hgm_waterchem/data_output/mdw_mod_GEE.rds")
unique(mdw.dat$PLOTNAME2) # 16 unique plot names but only 9 unique IDs
unique(mdw.dat$UCDID) # 16 unique plot names but only 9 unique IDs

# Violin Plots -------

# make a violin plot of NDVI
vio.ndvi<-ggplot(data=mdw.dat[mdw.dat$index=="NDVI" & mdw.dat$month>=6 & mdw.dat$month<9,], 
       aes(x=as.factor(month), y=mean, fill=as.factor(hgm_class_comb))) + 
  ggtitle("A Plot of NDVI") + xlab("Month") + guides(fill=guide_legend(title="HGM Type")) +
  geom_violin() + scale_fill_viridis(discrete = T) + 
  theme_bw() #+ facet_grid(month~.)

# make a violin plot of NDWI
vio.ndwi<-ggplot(data=mdw.dat[mdw.dat$index=="NDWI" & mdw.dat$month>=6 & mdw.dat$month<9,], 
       aes(x=as.factor(month), y=mean, fill=as.factor(hgm_class_comb))) + 
  ggtitle("A Plot of NDWI") + xlab("Month") + guides(fill=guide_legend(title="HGM Type")) +
  geom_violin() + scale_fill_viridis(discrete = T) + 
  theme_bw() #+ facet_grid(month~.)

vio_stack<-cowplot::plot_grid(vio.ndvi, vio.ndwi, align = "v",nrow = 2)
#cowplot::save_plot("./fig_output/violin_ndvi_ndwi_jun-aug.png", plot = vio_stack, ncol = 2,base_aspect_ratio = 0.7, base_height = 6)

violinplots

# Boxplots -------

# use notched boxplots to assess if medians are significantly different. Notches are used to compare groups; if the notches of two boxes do not overlap, this suggests that the medians are significantly different.

# make a box plot of NDVI
box.ndvi<-ggplot(data=mdw.dat[mdw.dat$index=="NDVI" & mdw.dat$month>=6 & mdw.dat$month<9,], 
       aes(x=as.factor(month), y=mean, fill=as.factor(hgm_class_comb))) + 
  ggtitle("Boxplot of NDVI") + xlab("Month") + guides(fill=guide_legend(title="HGM Type")) +
  geom_boxplot(notch = T) + scale_fill_viridis(discrete = T) + 
  theme_bw() #+ facet_grid(month~.)

# make a box plot of NDWI
box.ndwi<-ggplot(data=mdw.dat[mdw.dat$index=="NDWI" & mdw.dat$month>=6 & mdw.dat$month<9,], 
       aes(x=as.factor(month), y=mean, fill=as.factor(hgm_class_comb))) + 
  ggtitle("Boxplot of NDWI") + xlab("Month") + guides(fill=guide_legend(title="HGM Type")) +
  geom_boxplot(notch = T) + scale_fill_viridis(discrete = T) + 
  theme_bw() #+ facet_grid(month~.)

box_stack<-cowplot::plot_grid(box.ndvi, box.ndwi, align = "v",nrow = 2)

box_stack

cowplot::save_plot("./fig_output/box_ndvi_ndwi_jun-aug.png", plot = box_stack, ncol = 2,base_aspect_ratio = 0.7, base_height = 6)

boxplots

Use notched boxplots to assess if medians are significantly different. Notches are used to compare groups; if the notches of two boxes do not overlap, this suggests that the medians are significantly different.

# make a dotplot of NDVI
ggplot(data=mdw.dat[mdw.dat$index=="NDVI" & mdw.dat$month>4 & mdw.dat$month<9,], 
       aes(x=month, y=mean, fill=hgm_class_comb)) + 
  ggtitle("A Plot of NDVI") +
  geom_boxplot() + scale_fill_viridis(discrete = T) +
  theme_minimal()

# make a plot of NDWI
ggplot(data=mdw.dat[mdw.dat$index=="NDWI",], 
       aes(date, mean, group=hgm_class_comb, color=hgm_class_comb)) + 
  ggtitle("A Plot of NDWI") +
  # geom_text_repel(aes(label = hgm_classes), color = "black", 
  #                 box.padding = unit(0.5, 'lines')) + 
  geom_point() + scale_color_viridis(discrete = T) +
  theme_minimal()

Make a Shiny Map of Sites

library(leaflet)

m <- leaflet() %>% addTiles() %>% 
  #setView(lng = -120.8, lat = 39, zoom = 8) %>%  # set to Auburn/Colfax, zoom 5 for CA 
  addTiles(group = "OSM") %>%
  addProviderTiles("Stamen.TopOSMFeatures", group = "OSM Features") %>%
  addProviderTiles("Esri.WorldImagery", group = "ESRI Aerial") %>%
  addProviderTiles("Thunderforest.Landscape", group = "Topo") %>%
  hideGroup("OSM Features") %>% 
  

# proposed sites
  addCircleMarkers(data=mdw.dat, group="Weixelman Sites",
                   lng= ~LONG_DD, lat= ~LAT_DD,
                   popup=paste0("<b>Plot:</b> ", 
                                mdw.dat$PLOTNAME2, "<br>", "<b>HGM Type:</b> ",
                                mdw.dat$hgm_type, "<br>", "<b>UCD_ID:</b> ",
                                mdw.dat$UCDID, "<br>", "<b>Area_acres:</b> ",
                                mdw.dat$AREA_ACRE, "<br>", "<b>HUC12:</b> ",
                                mdw.dat$HUC12, "<br>", "<b>Elev_mean_m:</b> ",
                                mdw.dat$ELEV_MEAN),
                   stroke=TRUE, weight=0.6,radius=10,
                   fillOpacity = 0.8, color="black",
                   fillColor = "yellow") %>%
  
  # add controls for basemaps and data
  addLayersControl(
    baseGroups = c("OSM", "ESRI Aerial", "Topo"),
    overlayGroups = c("Weixelman Sites","UCD Mdws",
                      "OSM Features"),
    options = layersControlOptions(collapsed = T))

# Print Map
m